Prepare packages required for this analysis.

# Load Biobase.
if (! require(Biobase, quietly=TRUE)) {
  if (! exists("biocLite")) {
    source("https://bioconductor.org/biocLite.R")
  }
  biocLite("Biobase")
  library(Biobase)
}
# Load GEOquery.

if (! require(GEOquery, quietly=TRUE)) {
  if (! exists("biocLite")) {
    source("https://bioconductor.org/biocLite.R")
  }
  biocLite("GEOquery")
  library(GEOquery)
}

library(knitr)
library(edgeR)

Load normalized data from previous assignment.

normalized_count_data <- read.table(file=file.path("~/GSE125150_normalized_counts.txt"))
# Check if data imported correctly:
normalized_count_data[1:5,1:5]

Prepare for creating a heatmap and downstream analysis.

# Converting to matrix, because we have a data.frame:
heatmap_matrix <- normalized_count_data[,3:ncol(normalized_count_data)]
rownames(heatmap_matrix) <- normalized_count_data$gene_id
colnames(heatmap_matrix) <- colnames(normalized_count_data[,3:ncol(normalized_count_data)])
heatmap_matrix

library(ComplexHeatmap)
library(circlize)

# Create a heat map
if(min(heatmap_matrix) == 0){
    heatmap_col = colorRamp2(c( 0, max(heatmap_matrix)), c( "white", "red"))
  } else {
    heatmap_col = colorRamp2(c(min(heatmap_matrix), 0, max(heatmap_matrix)), c("blue", "white", "red"))
  }
current_heatmap <- Heatmap(as.matrix(heatmap_matrix),
                               show_row_dend = TRUE,
                               show_column_dend = TRUE, 
                               col=heatmap_col,
                               show_column_names = TRUE, 
                               show_row_names = FALSE,
                               show_heatmap_legend = TRUE
                               )


# Apply row-normalization
# Subtract the mean from each value and divide by the standard dev of the row to get row-normalization
# (scale() function does this for us).
heatmap_matrix <- t(scale(t(heatmap_matrix)))
if(min(heatmap_matrix) == 0){
    heatmap_col = colorRamp2(c( 0, max(heatmap_matrix)), c( "white", "red"))
  } else {
    heatmap_col = colorRamp2(c(min(heatmap_matrix), 0, max(heatmap_matrix)), c("blue", "white", "red"))
  }
current_heatmap <- Heatmap(as.matrix(heatmap_matrix),
                               show_row_dend = TRUE,
                               show_column_dend = TRUE, 
                               col=heatmap_col,
                               show_column_names = TRUE, 
                               show_row_names = FALSE,
                               show_heatmap_legend = TRUE
                               )

# Roughly, our data clusters based on treatment (KLF-1+ or KLF-1-), also shown by the MDS plot in A1.
plotMDS(heatmap_matrix, labels = rownames(samples), col = c("darkgreen", "blue")[factor(samples$Condition)])

Create a model, calculate p-values, and apply multiple hypothesis testing.

# Define groups:
samples <- data.frame(
        lapply(colnames(normalized_count_data)[3:10],
               FUN=function(x){
                 unlist(strsplit(x, split="_"))[c(2,1)]}))
colnames(samples) <- colnames(normalized_count_data)[3:10]
rownames(samples) <- c("sample", "treatment")
samples <- data.frame(t(samples))
samples

# Create data matrix
model_design <- model.matrix(~ samples$treatment)
kable(model_design, type="html")

expression_matrix <- as.matrix(normalized_count_data[,3:10])
rownames(expression_matrix) <- normalized_count_data$gene_id
colnames(expression_matrix) <- colnames(normalized_count_data)[3:10]
minimal_set <- ExpressionSet(assayData=expression_matrix)

# Fit the data to our linear model
fit <- lmFit(minimal_set, model_design)

# Apply empirical Bayes to compute differential expressions for the model
fit2 <- eBayes(fit,trend=TRUE)

topfit <- topTable(fit2, 
                   coef=ncol(model_design),
                   adjust.method = "BH",
                   number = nrow(expression_matrix))

# Merge gene names to topfit table
output_hits <- merge(normalized_count_data[,1:2],
                     topfit,
                     by.y=0,by.x=1,
                     all.y=TRUE)
# Sort by pvalue
output_hits <- output_hits[order(output_hits$P.Value),]
output_hits

# Check the results
kable(output_hits[1:10,],type="html")

length(which(output_hits$P.Value < 0.01))

## [1] 2841

length(which(output_hits$adj.P.Val < 0.01))

## [1] 1337

For Multiple Hypothesis Testing, I used the BH method. According to the documentation and literature, BH method controls for the false discovery rate (a less stringent condition than family-wise error rate) and therefore is a more powerful method than those that control for family-wise error rate. My dataset allows many genes to pass even with the stringent condition (with BH), so I decided to use BH. Before correction, 2841 genes passed a cutoff of 0.01, and after correction with the BH method, 1337 genes passed. With a cutoff of 0.05, there were 2859 (after correction) and I considered this as excessive, so I decided to make the cutoff more strict. I could have made the cutoff higher to 0.001; however, I was afraid I would disregard interesting genes, so I kept the cutoff at 0.01 as it is a generally accepted p-value for high significance.

Create heat map based on p-value cutoff.

top_hits <- output_hits$gene_id[output_hits$adj.P.Val<0.01]
top_hits
heatmap_matrix_tophits <- t(
  scale(t(heatmap_matrix[
    which(rownames(heatmap_matrix) %in% top_hits),])))
if(min(heatmap_matrix_tophits) == 0){
    heatmap_col = colorRamp2(c( 0, max(heatmap_matrix_tophits)), 
                             c( "white", "red"))
  } else {
    heatmap_col = colorRamp2(c(min(heatmap_matrix_tophits), 0, max(heatmap_matrix_tophits)), c("blue", "white", "red"))
  }
current_heatmap <- Heatmap(as.matrix(heatmap_matrix_tophits),
                           cluster_rows = TRUE,
                           cluster_columns = TRUE,
                               show_row_dend = TRUE,
                               show_column_dend = TRUE, 
                               col=heatmap_col,
                               show_column_names = TRUE, 
                               show_row_names = FALSE,
                               show_heatmap_legend = TRUE,
                               )
current_heatmap

Based on the heatmap, the conditions (KLF-1+ vs. KLF-1-) cluster together. All the samples in KLF-1+ (K2PT) share similar differential expression profiles and all the samples in KLF-1- (K2MT) share similar differential expression profiles.

Retrieve thresholded list for upregulated and downregulated genes

output_hits[,"rank"] <- -log(output_hits$adj.P.Val,base =10) * sign(output_hits$logFC)
output_hits <- output_hits[order(output_hits$rank),]

upregulated_genes <- output_hits$gene_name[
  which(output_hits$adj.P.Val < 0.01 
             & output_hits$logFC > 0)]

downregulated_genes <- output_hits$gene_name[
  which(output_hits$adj.P.Val < 0.01 
             & output_hits$logFC < 0)]

upregulated_genes[1:10]

##  [1] B4GALT6       MAP7D1        AKAP12        CLTCL1        ABLIM1       
##  [6] IL10          NDUFAF5       RP11-550P17.5 TRGV6         DNM1L        
## 13869 Levels: A1BG-AS1 A2M A2M-AS1 AAAS AACS AADAC AADAT AAED1 AAGAB ... ZZZ3

downregulated_genes[1:10]

##  [1] CCNA1   ARMC9   CD81    ALDH1A2 ID2     MAN1A1  PLP2    CTNNB1 
##  [9] H6PD    HS3ST1 
## 13869 Levels: A1BG-AS1 A2M A2M-AS1 AAAS AACS AADAC AADAT AAED1 AAGAB ... ZZZ3

write.table(x=upregulated_genes,
            file=file.path("KLF1_upregulated_genes.txt"),sep = "\t",
            row.names = FALSE,col.names = FALSE,quote = FALSE)

write.table(x=downregulated_genes,
            file=file.path("KLF1_downregulated_genes.txt"),sep = "\t",
            row.names = FALSE,col.names = FALSE,quote = FALSE)

ORA with g:Profiler

The extracted genes were put through g:Profiler for ORA. g:Profiler was used because of its easy manipulation. The following options were used. Significance threshold: Benjamini-Hochberg (Again, because of its powerfulness due to its stringence.) User threshold: 0.05 Annotation sources: GO molecular function, GO biological process, Reactome, WikiPathways (These sources were used because we want to know the mechanism of action behind the effect of KLF1+ macrophages on the maturation of erythrocytes, and what kinds of pathways may be involved for this maturation to occur. Other sources, like GO cellular component, were not relevant because we are interested in a broader picture on the pathways and processes)

Below are the images of related terms for the upregulated genes:

knitr::include_graphics("upregulated_genes_1.png")

knitr::include_graphics("upregulated_genes_2.png")

knitr::include_graphics("upregulated_genes_3.png")

knitr::include_graphics("upregulated_genes_4.png")

knitr::include_graphics("downregulated_genes_1.png")

knitr::include_graphics("downregulated_genes_2.png")

knitr::include_graphics("downregulated_genes_3.png")

Interpretations

The results of this analysis broadly supports the original paper’s conclusion, in that KLF-1 induction does lead to differentiation to macrophages that are involved in erythrocyte maturation. The result of g:Profiler shows that upregulated genes are involved in immune system processes, as well as cell-cell adhesion and communication. Also, protein kinase activity seems to be upregulated, and this makes sense because kinases are often related with signalling pathways and activations. Interestingly, these genes are also involved in neuronal processes which could potentially unveil interesting relationships between immune cells that are resident in umbilical cords vs. near neurons. Furthermore, the downregulated genes appear to be involved in kinase inhibition, which corresponds to the upregulated genes. One interesting process that is shown to be downregulated in KLF-1+ macrophages is the regulation of CD4 T cells. This could possibly be due to the fact that these macrophages are obligates to helping the maturation of erythrocytes.